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On Coulomb wave functions 


By Per O. M. Otsson 


Introduction 


The Appell functions [1], which are generalizations of the hypergeometric func- 
tions to functions of two variables, have been used in the quantum mechanical 
treatment of problems connected with the motion of charged particles in a Cou- 
lomb field. An extensive application is found in the study of nuclear structure 
by Coulomb excitation presented recently by K. Alder, A. Bohr, T. Huus, B. 
Mottelson and A. Winther [2]. 

From the fact that the matrix elements for electromagnetic transitions be- 
tween two states in an external Coulomb field can be expressed in terms of 
Appell functions in various ways and that these functions in many cases furnish 
analytic continuations of each other, their usefulness is obvious. These matrix 
elements may be looked upon as the Laplace transforms of radial Coulomb wave 
functions or products of such functions, and the many relations between different 
Appell functions allow a variety of expressions for the transition matrix elements. 
On the other hand the Laplace transforms of confluent Appell functions, when 
integrated over one of the variables only or with a simple relation between the 
two variables in order to form a function of one variable, often leads to expres- 
sions that are very similar to those of the transition matrix elements. This in- 
dicates a connection between radial Coulomb wave functions and some confluent 
Appell function. Since the Coulomb wave functions for non-bound states are 
rather cumbersome to work with it was considered worth while to look for a 
representation of the radial Coulomb wave functions in terms of a confluent 
Appell function. There were reasons to believe that the latter function, in some 
respects at least, should have more attractive features than the Whittaker functions 
ordinarily used in this connection. An investigation based on these arguments is 
presented below. 

Real valued wave functions have been chosen throughout. This is not a ne- 
cessity, but the symmetry gained thereby seems to make it desirable. 

Although the transition matrix elements between two Coulomb states are con- 
veniently expressible in Appell functions there remains the difficulty that the 
power series expansions of the Appell functions, which perhaps provide the simplest 
way of obtaining numerical values except possibly for very special cases, often do 
not converge or converge very poorly for practical purposes. 

The functions introduced in the present paper allow more transformations than 
the ordinary Coulomb wave functions to which they are reduced for special values 
of the parameters and seem to offer more possibilities to extend the domain of 

mvergence. 
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Radial Coulomb wave functions 


We have chosen to treat the relativistic case but there is no difficulty in 
translating the derived formulas to their non-relativistic analogies. Since the an- 
gular dependence of a partial wave is the same for all potentials with spherical 
symmetry and this angular dependence may be found elsewhere [3] we shall be 
concerned here with the radial functions only. They satisfy the differential equa- 
tions [3] 


df(r) x-1 ey ts a Z 
smears, (z Ite 2) a(n) 


d 

d - Z +1 ) 
ANA SS eae! Rad Se 

ip - (E414 2) F(0 E g(r), 


where E is the energy measured in units of mc”, r is measured in units of h/me 
and x is connected with the total angular momentum number by the relation 
+u=j+1/2. 

It is convenient to introduce new variables 


4 (r)=rf (r) 


and sometimes also a new independent variable t=2«Zr, which leads to the 


equations 
Ei— 1 Z 
dx, %y a: ( a Lo = 0 


diakinl oo fae 
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We claim that the equations (3) have solutions of the form 


aZ E-1 d® (y, EL, i] 
= if | ——— = ee Se eee 
%, (t) oot 9 ® (y, E, t) + (y—%) di d 
E+1 dO (y, E,t 
X(t) = -r| O (y, H, t)+(y +x) mE), 
where y= 4+ Vie — eB 
and ® (y, H,t) is any solution of 
@O(y,E,t) 2ydO(y,H,t) (EH a . j 
dt? Qa arr: P aca Die 6) 


This is easily verified if the expressions (4) are introduced into the equations 
(3). Both expressions on the left hand side of (3) are then found to be identi- 
cally zero if ® (y, Z,t) satisfies (5). This equation has, independently of the sign 
of y, always a solution that is finite for t=0 and which is uniquely determined 
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by its value at this point. We shall denote by ®(y, #,t) the solution of (5) for 
which 


® (y, H,0)=1. 


It follows that two independent solutions of the equations (3) are obtained from 
(4), corresponding to the two possible signs of y, and the problem of solving (3) 
is reduced to finding the regular solution ® (y, H,t) of (5). 

The solutions @(y,H#,t) of the equation (5) are well known and are expo- 
nential functions times confluent hypergeometric functions. We wish, however, 
to point out that the regular solution of (5) may be conveniently represented 
by a confluent Appell function [1]. The function in question is 


48,8: 27; 2 y)— 5 Pim Bazi y” | (6) 


The double series is convergent for all values of x and y independently of the 
values of the parameters with the exception of 2 =negative integer or zero. This 
exception is obvious from the definition of the symbols (f)m etc. which are ab- 
breviations for products of the type 


(B)m =B(B +1) (B +2)... (B-+m—1). 


In order to obtain a function of one variable we put x=qt and y=st: 


©, (6, 6's 2y; at, et) = $ (B)m(B’n g™ 5" 


k=0 (2), mfn=k m!n! 


‘The second sum may be expressed in a hypergeometric function: For [4] 
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Also more directly 
(l1—gt) 8 (1—st)-" = S y (B)m (B'n g” 8” 


k=0 m+n=k m!n! 


(B)n(B')n 8" _(B+B’e a 


min=k m!\n! k} 


Hence 


oF, (—kp' 8+ 61-2) 


and consequently 
f 


er SBT BW pifac: gi moih_® hs 
©, (8,8';2y;qt, st)= > (2y)n 2! £F,( n;B ;B+B 31 *) (7) 


To get a solution of (5) we make the following choice for the parameters 
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where n= ,p=VE*-1 


and obtain 


S CPT) FB (—n,y —ins 2932) (9) 


n=0 n! 
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This function is real valued and solves equation (5). The first statement follows 
directly from the definition of the function and the second statement can be 


proved as follows. 


Put O(y, H,t)= de : 


, nl (10) 
and use equation (5) to determine the coefficients c,. The following recursion 


formula will be obtained 
2 


(n+ 2y)eni1+ Bent pia gaten _1=0, 


which is satisfied by the real coefficients 


ptt ee Oy nn clare (11) 
n 20 Z, med +Y q; Y5 
since from the recursion formulas for the hypergeometric functions it can be 
shown that 


(n+2y).Fy(—n—-ly-tns2y; 2)—2ingF,(—n, y—195 273 2)— 
—nF,(—n+1,y—in32y;3 2) =0. 


Ne © (y, Bt) 5 OPO Fy (=m, yi 2752) 


which is the function (9) and thus 
® (y, E, t) =O, (y +in,y— 19s 2st pr, —tPr) (12) 


The series expansion (10) for ® (y, E,t) is convergent for all values of t and # 
as a consequence of the afore-mentioned convergence of the @,-function of 
which it is a special case. The coefficients ¢c, are polynomials in #, as can be 
seen at once from their recursion formula, and therefore ® (y, H,t) is an ana- 
lytic function of both the variables E and t for any finite value of these variables. 
There is no need then to distinguish between bound and non-bound states from 
a formal point of view. 
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The function (12) is real and ip enters only in even powers. As mentioned 
before it supplies us with two independent solutions of the Dirac equation through 
(4). The solution that is regular at the origin is obtained for y= +Vx?—o2Z? 
and the irregular solution for y= —Vx?— 0? Z?. 

It has been mentioned already that ® (y, #,f) must be an exponential function 
times a confluent hypergeometric function and it is not difficult to show that 


O (y, E, t)=e'?" Fy (y+in;2y;2ipr), (13) 


e.g. by solving (5). We will, however, prove (13) in a different manner and at 
the same time arrive at a more general result, that will be useful in applications. 
The notation pF, (a1, %,...%p3 Py, Bs,... Bg; x) will be used throughout this paper 
for a generalized hypergeometric function. They are defined by the power series 
expansion 


pg (ot, a, «++ 5 By, By, -.. Bai &) eal = eae nee ai 


and its analytic continuations. 
The proof of (13) is based on the Appell function [1] 


- (2) min (B)m (B’)n Oe 


F, (a, B, B’; 2; 2, y)= (2Y)minm! n! 


(14) 


n=0 


from which the confluent Appell function ®, is obtained by the limiting process 
of confluence: 


, : 1 ; 
2(B, P';2y;«,y)=lim Fy (7.8.8 2:62.64) (15) 
: 1 
as lim (*) (ex)™(ey)*=a7y"- 
e>0 \E/ min 


Before computing the limit in (15) a well known transformation of the F,- 
function is applied [1]: 


f -@ ’ ar x mm 
F, («, B, B’;2y;2,y)=(l—2) * Fy («.27-6-'.p 32737 y="). 


x 
Then 
®, (8, 8’; 232, y) = 
Sra iasan\ =i" _é% fet) 
So - A; Ff de Oe —1 l-ex 
=e", (2y —B—f’; B’; 2y; —x,y—2) (16) 


4 
and if we put fb’ +B=2y in agreement with (8) the ®,-function on the right 
nd side is reduced to a confluent hypergeometric function 
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®, (0, B';2y;2,y—x) =F, (B53 2y;y—2) 


This is a direct consequence of the definition (6) of the ®,-function and proves 
that the @,-function is reduced to an exponential function times a confluent 
hypergeometric function if 6+)’ =2y or 


®, (8, 8’; B+ P's x, y) =e" F, (B'; B+ B's y—2). (17) 
Also for any 4 
e’ , (B, B’; B+ B's x, y) =D_(B, B's B+B';xt+Asy+A). 
Using the notations (8) we obtain 


O,(yt+in,y—in;2y;tpr,—ipr) =e" FP, (y—in;2y;—2ipr)= 
=e '?' Fi (y+in;2y;2ipr), 
where the last equality follows from the fact that the function is real valued. 
This proves (13) and connects our formalism with more familiar expressions for 


the Coulomb wave functions. These functions are usually given in terms of the 
function 


H (r)=e'?" LF, (y+1+in; 2y+1;2ipr) 


This function may be expressed in terms of the functions ® (y, #,t) with the 
aid of recursion relations for the ,/,-functions: 


ReH=}(H+H"*)=3e '""(,F,(lt+ytin;l+2y;2ipr)+ 
+,F, (yt+in;l+2y;2ipr]=pe'?" | Fiybinsd +2ys2ipn + 


2ipr 
1+2y? 


PL by + ins2y+2:2ipr) + Fy Gy bins L427 2ipr] 


, P 2ipr . 9:93 
eet: .F: (bin b2y:2ipe) tye Fey indy +2i2ipr| 


Et : a 
=e?" [uly +ins2ys2ipe) + oh ty +inidy +2;2ipn) 
2 


Et 
—— O(y +1, EZ, t), (18) 
(2y)o 7 


and for the imaginary part after similar calculations 


=O (y, H,t)+ 


Im B= 5G) Py +1 Bt. (19) 


Since generally 
d, F, (a; ¢; x) 


a 
de =-1F (a+ le+1;2) 
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it follows that 


dt  2aZ 2aZ 2y 
and since On 80 is real valued 
d® (y, E, t) p EH 
i sg mi ieee 


and with the aid of (18) and (19) finally 


d® (y, EL, t) 


4 rey 


=E®(y, B,t)+ © (y +1, E£,t). (20) 


(2). 


It is sometimes an advantage not to have the derivative of O(y, H,t) in the 
wave functions. With the aid of (20) it can be eliminated yielding 


aZ Hu-y Livesey t 


uty (2y)s 


“-y 2y 


1 


t’ Jo (y, EB, t)+ M(y+1, #, | 


Eut+y ¢ 
WAS) )s 


L,= ote log B+ D(y+1,8.0| 


For the sake of completeness we will also derive the normalization constant 
for the wave functions for non-bound states. The normalization constant for 
bound states is found e.g. in reference [3], and correcting it to fit the solutions 
(21) is a very simple operation. The solutions are normalized to a flux of one 
particle per unit energy interval per second [5] 


z— 
Nh / = 


= sin prin In 2pr—2L— arg (y bin) +0| 
E+1 

N, al, 

the ap 


where 6 is an asymptotic phase. This asymtotic form is most easily obtained 
from the well known asymtotic behavior of the function ,F, (y+in;2y;2ipr) 
for large arguments. 


(22) 


cos |pren In 2pr— "2 argT (y-+in) +9} 


; : ety em eztpr 
1 (y+in;2y;2ipr)~T (2 i Formentera T (y+in)(2ipry al 
Using (13) we obtain 


2D (Qype 2 
IT (y+in)| Qipry 


a 


| O(y, £, t) ~ cos | prt In 2pr— 22 arg (y-+in)| 
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aD(y,H,t) pp T'(2y) es 
dt a&|V(y+in)| (2ipr) 


> Sin aw In 2pr— arg Ty +in)| : 


These expressions combined with (4) give the asymptotic form (22) where 


en 218 ytin 
iaZ 
ee 
Pp 
Z(E-1) 
or ta 6 ele 93 
¢ p(%—y) eas 
Gay [Potinle® 
p ‘= IT (y tin) le? | 
d N,= D4 
ne i (5) Ex-y V2x pT (2y) at 


The second expression for the asymptotic phase is a consequence of the first 
formula. In the real valued formalism presented here the first expression never 
will be needed but it is the definition of 6 usually encountered in papers dealing 
with Coulomb wave functions. Writing down the normalized wave functions we have 


Hut+y t 
= Y 
NX, =aZNt J (B04 Se ay oy th Bo] 
Ex+y t 
Ny Xo = (*#— velo , H, t)4 ® +1,8,0]| 25 
Xe = (x—) (y, £, t) er ONS (y ) (25) 
eae | an 
p ‘ye? (y +1) |e? 

h = : 
ee 2 BB u—y V2apl (2y+1) 


In deriving (25) care has to be taken not to loose a change of sign which 
shall be there when the quantum number x changes its sign. 
The quantity under the square root sign is positive for all occurring values of 


the three parameters involved. It replaces in a certain sense the factors V1-E 


and 1+ in the ordinary expression for the relativistic Coulomb wave functions, 
but here the change from bound to non-bound states introduces no imaginary 
quantities. It should perhaps be remarked that the normalized irregular wave 
functions may still be obtained by changing the sign of y wherever it occurs. 


Integral representations 


Integral representations for the ®,-function valid for any values of the para- 
meters are possible only as contour integrals in the complex plane. Such a re- 
presentation has been given by Erdélyi [6] and a representation in the form of 
a double integral over a real domain by Appell and Kampé de Feériet [1]. Inte- 
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gral representations are useful tools in the transformation theory of the function 


represented and also in connection with analytic continuations. 
We will show that the ®,-function has an integral representation 


®, (8, B’;2y; x,y) = 


1 
P(2y) | p-1 2y-B- 
tut (1 —u)y?? oF [B's 2y—Bsy(l—u)] du (26 
The integral has a meaning only when Re2y>Ref>O0 but results derived 
from it are in most cases free from this restriction. To prove (26) we write 


1 
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=0(2y—B),n! T(2y+n) Fy (Bp; 2y +n; 2) 
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[P(2y+n) zu, pl Se Voge te p=1 
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This is a well known integral representation of the ,F,-function. Then 
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pel (y) zu, p-1 __ ,\2y—-B-1 Ge ah ys ze = 
TR a I et i nce 


2 (B ny” erie 
= cea is (B;2y+n;2)= 


_ S (Bn (B)my” =” ge} 
e aviascminl O2(8, Bs 252, 9). 


For the choice B+f/’=2y we obtain 
1 


; D (BBB+ Ba N= Rare [errr tau tae (28) 
0 
‘since Fy (6'; B'sy (1-w) =e". 
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This shows, as we already know, that the ®,-function with the indicated relation 
between the parameters is of the type 


®, (6, B's B +B’; x, y) =e" f (w—y) 
Using (27) we obtain 


®, 6, B’ B+’; x,y) =e" 1F,(B;B +p’; x—y) (29) 


and for the parameters (8) 


—ipr 


DO, (y+in, y—in;2y;ipr,—tpr) =e Fy (y+in;2y;2ipr) 
in agreement with the previous result. The interesting result also previously 
derived 


®, (8, B's B+’; 2, y) =e 7 ®, (B, B's BR +B sx+A,y+A) (30) 


is an easy consequence of the integral representation and is an example of an 
useful transformation, in that it allows a certain arbitrariness in the choice of 
the arguments of the function which is of importance in integrals over Coulomb 
wave functions where this freedom may be used to select convergent expansions. 


Expansions of the wave functions 


Besides the power series expansion, the general term of which is derived in (11), 
a number of expansions of the wave functions in terms of hypergeometric 
functions may be found. The term hypergeometric function is here to be 
interpreted in its more general sense. Expansions of this kind may be obtained 
by expanding the Laplace transform of the wave function and then apply the 
inverse Laplace transform to the terms. 

We have, integrating term by term, 


L{P1* @,(B,B';2y; xt, yt)} = [ e*'P1 1, (8, Bs 2y; xt, yt) dt= 
0 


See (1 -2)" (1 -)" (31) 
s?? 8 8 


2 a\ ~8 y = ih¢ °° l 
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where z is an arbitrary parameter, the Laplace transform can be written 


Expanding 


z 
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Li?’ @, (B, B's 2y3 et, yb} =D (2y) 3 dn ayaa (32) 
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where the definition of the )/,-function is as in (13). This function is essentially 
a Bessel function. The relation is 


J, (Z)= Tae Ie Ae 


T(v+1)° 3 
From (32) and (33) it follows that 


L{t?’ 1 ®, (B, 8’; 2y; xt, yt)}}= Deas Prind F Qyt+n; -20)| 


(2y)n 


from which we conclude 


OatB Pavia y= $s 


” F,(2y +n; —2t) (34) 


when the right hand side is uniformly convergent. 
From (32) we obtain 


Ba (B)x (Br x y'2™ 
as k+l+m=n kil! m! 
which we shall prove to be equal to 
(B')ny” ( Cdn ‘| 
na —N,;P; Ihe Tae tee 35 
a ay arn, 6 fp —7 #5 (35) 
where the ©,-function is one of the confluent Appell functions. It is defined by 
bse: fGen v | 36 
®, (a, Ont, ) 2 ape! ( ) 


Using this definition in (35) we obtain 
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Then 


pe Seepead idl ys See 
0, (=n, 61-8 Be 5) ee? k! 1! m! 


and (35) is proved. 

An interesting feature of the expansion (34) is that the dependence of the 
energy H can be separated from the function depending on r and thus be con- 
tained in the coefficients only. In problems where a knowledge of the wave 
functions in the neighbourhood of the origin is needed for not too high energies 
the expansion should be quite useful. 

For 6+ f'=2y the left hand side of (34) is an arbitrary exponential function 


times a radial Coulomb wave function, and if also «= -y=5-5 the function 
® (y, £,t) is obtained: pe 


Ply. Hf dita, yt yon a ites ce (37) 
(5?) 
where an = nae bee aN (38) 
2aZ 


Expansions of confluent hypergeometric functions have been given by F. Tricomi 
[7] and one of his expansions is obtained from the one given here (37) for z= £. 
Tricomi gives a recursion formula for the coefficients in this case which in our 
notation is 


2 
(n+1) anys = —piegallnt2y— lays — Baya] (39) 
2 
P 
a) = 1, a,=0, a= aR 
Ep (y)2p4 


etc. 


73 19927? “4 39 A Ze 


The coefficient a, is a polynomial in p with possible terms 


Ono, Dib tgon Pet +en-pajp" ls] (40) 


where 5 | means the largest integer < s but no odd powers occur. As an example, : 


the coefficient a,; is of the form 
Ay5 = Cy PY? + Cy p™ + Cyq pt 
according to (40) but with all odd powers of p omitted. 
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The series (34) contains (7) as a special case for z=0 since then 


ay = On", (—n, Bs 1-32, 0) eal | 


Using a well known transformation of the ,/,-function we obtain 


n ae a) 
oy DFP, (= 


n, Bs B+B'; 1-2) 


which gives the expansion (7). The general expansion (34) thus contains the 
power series expansion as well as different expansions in cylindrical functions and 
has been derived mainly in order to study the energy dependence of matrix 
elements for transitions between Coulomb states. 

A direct use of the Tricomi expansion, e.g. in the matrix elements for electro- 
magnetic transitions from a bound to a non-bound state, leads to an expansion 
that converges only for very small photon energies. The more general expansion 
given here and particularly the use of the Appell functions greatly increases the 
number of representations of the Coulomb matrix elements. This may lead to 
more convenient ways of computing them. As an example the radiative transition 
matrix element from a bound to a non-bound state as occurring in internal conver- 
sion can be expressed by means of the Appell function F, and the series repre- 
sentation of this function will then converge for all photon energies. The double 
series expansion can be expressed as a single series. For high photon energies the 
convergence is slow but tentative numerical investigations indicate that the ex- 
pansion is practically useful for photon energies as high as ten times the binding 
energy of the emitted electron. 
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Communicated 12 November 1958 by Oskar KtLEtrn and Ivar WALLER 


The diffusion cooling effect in heavy water 


By N. G. Ssy6stTRAND 


An attempt has been made to determine the diffusion parameters of heavy water 
using the pulsed neutron source method [1]. Because of the rather low counting 
rates no complete investigation was possible. However, some results were obtained 
for the diffusion cooling coefficient and they are reported here. 

Two quantities of heavy water were used: 3.68 and 3.14 kg with an H,O impu- 
rity of 0.60 +0.10 and 0.75 +0.10 per cent by weight respectively. The heavy water 
was contained in a cylinder 18 cm in diameter, and the experimental arrangement 
was essentially the same as in earlier work [2]. The temperature of the moderator 
was 20°C. Measurements were made with different pulse lengths and channel widths 
and with different detector positions. 

The bucklings for the two geometries were calculated as 0.0864 and 0.0966 cm-?. 
The uncertainty is about 0.002 em-, and it arises from the estimation of the extra- 
polation length. For pure D,O at 20°C the transport mean free path was assumed 
to be 2.51+0.04 cm [3,4,5]. The value adopted for the extrapolation length was 
1.7+0.1 em, taking the impurity and decreased neutron temperature into account. 

The measured time constants for the two bucklings were 14640 +200 and 16 840 
+400 sec-!. Subtracting the small contribution from absorption, dividing by the 
bucklings and correcting for the H,O impurity, the diffusion constants for pure 
D,O are 1.74+0.06 and 1.80 +0.08, both in units of 105 cm?/sec. The value related 
to the smaller buckling is smaller than the other, but considering the limits of error 
this is not unreasonable. 

In order to calculate the diffusion cooling coefficient c we use the relation 


D=D,-— (c—d) B?, 


where D, is the diffusion constant at zero buckling and d is a correction term [6] 
which for heavy water is estimated to 1.3x 10° cm/sec. Using a value of (2.08 + 
0.03) 10° em?/sec for D, [3, 4,5], c—d is found to be (3.5 +0.8) 10° cm4/sec, and the 
diffusion cooling coefficient is (4.8 +1.0) 10° em4/sec. This may be compared with 
the theoretical value 8.3 10° cm4/sec from the assumption of a monatomic gas [7] 
with nuclei of mass 20. The discrepancy indicates that the average energy transfer 
between neutrons and moderator is somewhat larger than the theory predicts. 
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Communicated 12 November 1958 by Oskar Kiern and IvAak WaLLEer 


On the theory underlying diffusion measurements with 
pulsed neutron sources 


By N.G. Sy6sTRAND 


SUMMARY 


Starting from the time-dependent Boltzmann equation some problems relating to measure- 
ments with pulsed neutron sources in small geometries are investigated. It is shown that the 
commonly used formula for the time constant of the fundamental mode of the neutron distri- 
bution is only approximate. As a consequence all diffusion cooling coefficients measured up to 
now are too small by 5 to 30 per cent, Further, the conditions under which the results from 
pulsed neutron source measurements in small geometries can be applied to a reactor are given, 
and the boundary conditions to use in the interpretation of such experiments are discussed. 


Introduction 


During the last few years many investigations of the neutron diffusion proper- 
ties of moderators have been made with the pulsed neutron source method in 
a small geometry [1-10]. The time decay constant A, of the fundamental mode 
of the neutron distribution is then given by 


Ao=Aa+ DB, (1) 


where /, is the absorption probability per unit time, D the diffusion constant for 
neutron density and Bj the geometric buckling. As the neutron temperature va- 
ties with buckling due to the diffusion cooling effect and D is temperature de- 
pendent, we can write in the first approximation 

D=D,-— CB, (2) 
where C is the diffusion cooling coefficient. 

_ In the experiments the time constant J, is measured as a function of the buck- 
ling. By fitting the data to the parabola; 


Ip = Aq + Dy BY — CBS (3) 


it is possible to determine the parameters 1,, D) and C. (For details of the method 
see reference 11). Conversely, if the paradheters are known for a moderator, eq. 
(3) can be used for determination of the geometric buckling. Such measurements 
n give useful information for reactors of peculiar shapes, on the influence of 
ntrol rods etc. [6, 9]. 


2 147 


N. G. SJOSTRAND, Diffusion measurements with pulsed neutron sources 


Eq. (1) is based on the ordinary diffusion approximation with all its limita- 
tions. In spite of the large experimental effort in this field, no thorough discus- 
sion of the theory underlying the pulsed source experiments has appeared in the 
literature before. In his book on neutron transport theory Davison [12] shows 
that the time-dependent Boltzmann equation can be reduced to the same form 
as the stationary equation, and he thereafter only treats the stationary case in 
detail. Even though many of the results obtained for the stationary problems 
can thus be applied directly to the time-dependent ones, the points of interest 
are so different in the two cases that it is worth-while to go into some detail. 

The questions to be considered are the following: 

1. Under which conditions can results obtained with pulsed source measure- 
ments in a small geometry be applied to a reactor ? 

2. When is equation (1) valid? 

3. What are the correct boundary conditions to use? 

Reference will constantly be made to the work of Davison [12] by the abbrevia- 
tion D. followed by the section number in his book. 


The equivalence between the stationary and time-dependent problems 


As the discussion following will be mainly based on the equivalence between 
the time-dependent and stationary problems, we first recapitulate the basis of 
this equivalence (D. 3.2). 

We consider a single, homogeneous, finite body, containing no sources, and 
assume also that no fissionable material is present. The Boltzmann equation then 
takes the form 
éN (r, vQ, t) 


apes + vQ grad N (r, vQ, t) + v Dtot (v) N (r, vQ, #) 


= | v' d, (v’) dv’ [fa (r, v' Q’, thf (v’ QD’ > eQ) dQ’. (4) 


Here N (r, vQ, t)dV dvdQ dt is the number of neutrons at time ¢ in the volume 
element dV around the point r, travelling with speed between v and v+dv ina 
direction lying within a solid angle d& around the direction 8. With our assump- 
tion of no fissionable material the total cross section Trot is given by Dror = Us + Le. 
where &, is the scattering and &, the capture cross section. The scattering function 
f (v’ QD’ +vQ) must be normalized so that 


[Pf (@ 2! > 2) dvdQ=1. (5) 
Equation (4) is valid for the neutrons in a moderator after the injection of 


burst of neutrons. If we make a Laplace _transformation of this equation with 
respect to time, we get an equation for N (r, v&,/), where 


N (r, oQ, a= [ N(r, vQ, the” dt. (6 

0 
The resulting equation is identical with the stationary form of eq. (4), if th 
macroscopic capture cross section is increased by A/v, and of course N (r, 0Q, t 


is replaced by N (r, vQ, A). 
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As there is no source term in the equation, we have an eigenvalue problem, 
because N(r,vQ,/) has to satisfy certain boundary conditions. Our assumption 
of a finite body makes the eigenvalues form a discrete set. The largest eigenvalue 
A, is the time constant of the system, and as we have assumed that no fissile 
material is present, it is necessarily negative. It follows that the decay of the 
neutron population sufficiently long after injection is exponential with the time 
constant Ap». 

Once we know this, it is easy to verify the above result on the apparent in- 
crease (or rather decrease, as A, is negative) of the capture cross section by the 
substitution 


N (r, vQ, t) =N (r, vQ) e* (7) 


in eq. (4). 

In the experiments which are of interest here, fast neutrons are injected into 
a finite moderator. The neutrons slow down to thermal energies and reach an 
equilibrium state, characterized by a constant spatial and velocity distribution, 
apart from the exponential decrease of the amplitude with the time constant Ap. 
It may be possible to find systems where this condition is not exactly fulfilled 
during times suitable for measurements. However, for most moderators the slowing 
down times are short compared to the life times of the thermal neutrons, and 
it is thus possible to await the thermal equilibrium and also to eliminate the 
higher harmonic modes. In fact, a pure exponential decay has sometimes been 
found experimentally over four decades. In the following we shall therefore only 
be concerned with the fundamental mode and its time constant A). 


The leakage and absorption rates 


If the fictitious capture is included in the stationary Boltzmann equation this 
can be written 


~Q grad N (r, vQ)+[—A,+v (2. (v) +B; (v))] N (r, vQ) 
=f v'X, (v') dv! [f N (r, v' Q’) f (v' Q'>vQ) dQ’, (8) 


where we have changed the sign of A, to get it positive. Integration over all neu- 
tron directions Q, all velocities v, and over the whole volume V of the body gives 


USE awed oS dodray [fff] oo grad Wir, 0) dodn.av 


os We. Q)dvdQaV =O] Nir oQ)dvdQav 


This follows because the integral of the right side of eq. (8) over v and Q is 
equal to 
: [ff eX. (v) N (r, vQ) dvdQ 


cording to eq. (5). 

The first term in eq. (9) is the absorption probability per unit time, which we 
ve denoted by 2, in eq. (1). If the absorption cross section &, follows the 1/v 
w, then 4,=v-,(v) independent of the neutron spectrum. This fundamental 
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result has been used for cross section measurements [13]. If decay curves are 
measured with and without an extra absorber and the results are divided point 
by point, one gets an exponential with the time constant v-¥{(v) where X¢ is 
the capture cross section of the extra absorber. As the spatial neutron distribu- 
tion generally is very little changed by the addition of an absorber this expo- 
nential is obtained also in the presence of higher harmonic modes. Thus apart 
from very small corrections for the change in the spatial distribution of the neu- 
trons due to the added absorber and for finite resolution in the measurement, 
the exponential behaviour is quite general even for rather short times after the 
neutron burst. 

If the capture cross section does not follow the 1/v law the first term in eq. 
(9) will depend on the neutron spectrum. This fact has been used to study the 
effective neutron temperature as a function of the buckling of the geometry [14]. 

In the second term of eq. (9) vQ grad N=div (vQN). Then Gauss’ theorem 
shows that the numerator is the net outflow of neutrons through the surface of 
the body. This term represents thus the leakage probability per unit time, and 
corresponds to the second term on the right side of eq. (1). 


The equivalence between a reactor and a pulsed moderator 


We assume now and in the following that the capture cross section follows 
the 1/v law. The time constant is then 


Ap=VUcty, (10) 


where » is the leakage rate. We introduce the quantities 


1 
ah case Pree (11) 
bins as) gn 
c*(v)= Se (12) 


and can then write equation (8) as follows. 


vQ grad N (r, v8) + ee 


= } mee i N(r, v' Q')f (v' Q' + rQ) dQ’ (13 
*(v') 


This equation is identical with that for the neutron density in a critical reacto1 
if instead of the quantities /*(v) and c*(v) we use the ordinary total mean fre¢ 
path and number of secondaries per collision respectively. They are defined ai 


1 


= Se) Be EH(0) a 
_ __ &s(v) +9 By (0) 
OW) =>, ih Sale aeien Ue 
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where &, is the fission cross section and v the number of neutrons emitted per fis- 
sion. Also the function f(v’ 2’—vQ) has a different form in the presence of fission 
neutrons (D. 2.2). Eq. (13) can be considered as an eigenvalue problem in c(v) 
or c*(v) and we see that both for a critical reactor and a pulsed moderator c must 
be larger than unity. 

The important conclusions from the above equivalence can be stated as follows: 

We have a bare, homogenous, critical reactor characterized by c(v),l(v) and 
a scattering function and a pulsed moderator characterized by c*(v), 1*(v) and a 
scattering function. The eigenvalues c(v) and c*(v) are identical, if 


(a) the functions f(v’ Q’—vQ) in the two systems are the same, 
(v)/l* (v) = const., 
(c) the ratio of the linear dimensions of the two systems is equal to this ratio 
between the total mean free paths. 


This result, though very general, is of little practical importance, as conditions 
(a) and (b) are seldom fulfilled. We will therefore try to find similar conditions 
for the one-group approximation (called by Davison the constant cross section 
approximation). 

Two assumptions are made. First, the neutron spectra in the systems are sup- 
posed to be independent of position. We can then write 


N (r, vQ) =p (r, Q) F (ev). (16) 


Secondly, the angular distribution of neutrons resulting from collisions is supposed 
to be independent of the velocity of the neutrons. This means 


f (v Q! > VQ) = g (v' > v) f (Q’ > Q). (17) 


Making these substitutions in eq. (13) and integrating over all v, we obtain 
(D. 4.4) 


Q grad y (r, Q) [ere v)dv+yp(r, Q) om 
=| [verano > Q) dQ’ ee (18) 


where we have omitted the asterisks. We may define averages / and c of | (v) and 
c(v) by means of 


1 fi nee Fame 
i | ered [AO : (19) 
0 0 
Pv F(v)dv [ae ) F (vw) dv Ye 
ef iw) J Uv) 
; 0 0 
Thus eq. (18) is finally reduced to 
: Q grad p(r, 9) +2" ETL wr, Q') f (Q’>Q) dQ’. (21) 


N. G. SJOSTRAND, Diffusion measurements with pulsed neutron sources 


This equation and the definitions of the averages are the same as Davison 
found but are more generally valid as we have not assumed the scattering to be 
isotropic. 

If the definitions (11) and (12) are used in eqs. (19) and (20) we can express 
l* and c* in terms of the average scattering cross section. We choose to normalize 
the neutron density spectrum, that is 


J F(v)dv=1. (22) 


2-3 ( v F (v) d, (v) dv (23) 
iy 
il 
we find then Paso, (24) 
x.— y/o 
one (25) 
x;— y/o 


where @ is the average neutron velocity. For a critical reactor the expressions (14) 
and (15) are unchanged when the averages are introduced instead of the velocity 
dependent quantities. 

If cis close to unity the effect of anisotropic scattering on the solution of eq. 
(21) is in first approximation equivalent to changing the total mean free path 
to the transport mean free path (D. 17.3). We can then formulate the conditions 
for getting the same eigenvalues ¢ and c*, namely 


(a) the neutron spectra in the two systems must be independent of position, 

(6) the angular distribution of neutrons resulting from collisions must be inde- 
pendent of velocity, . 

(c) the eigenvalues c and c* must be close to unity, 

(d) the linear dimensions of the systems must be the same when measured in 
units of the respective transport mean free paths. 


The equivalence between diffusion length and buckling 


We are now going to study the solution of equation (21), and we start with 
assuming that the scattering is isotropic. It has been shown (D. 5.1) that for an 
infinite medium a solution to the corresponding integral equation can then be 
obtained if 
L L+l 
= [lore 26 
CMG °8 L-l gs 


and the solution satisfies the equation 
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V?e(r)=9(r)/L, (27) 
where o(r)=|{ pir, Q)dQ. (28) 


The diffusion length Z is a real quantity if c<1. If c>1, LZ becomes imaginary, 
but the derivation given by Davison is obviously valid also in this case. Introducing 


1 
ae 2 
B TZ (29) 
A : Tye NL oF 
we find the relation eee y tan Bl (30) 


and the solution of the Boltzmann equation satisfies then the equation 
V?0(r) + Bo (r)=0. (31) 


Even if we have derived these formulae with the assumption of an infinite 
medium we can also apply them approximately to finite systems if suitable boun- 
dary conditions are imposed. However, the flux o(r) calculated from this diffusion 
approximation will as usual be decidedly inaccurate within a boundary layer se- 
veral mean free paths in thickness. 

Obviously the above relation between B, J and c¢ (eq. 30) is equally valid for 
both physical cases for which ¢>1, namely for a critical reactor and for a pulsed 
moderator, if only we use the proper values for / and c. In both cases eq. (31) 
gives the neutron distribution, and the lowest eigenvalue determined by the bound- 
ary conditions is called the geometric buckling, B?. 

Eq. (30) can be interpreted for a reactor as giving the material buckling By, 
for a medium with the properties defined by ¢ and 1. For a critical reactor Bj, = 
B}. This condition is already inherent in eqs. (30) and (31) and the boundary condi- 
tions because eq. (21) describes a source free, stationary state, i.e. criticality. 


The relation between time constant and buckling for isotropic scattering 


From eq. (30) we can now find a relation between the time constant and geo- 
metric buckling for a small moderator. By substituting /* and c* from eqs. (24) 
and (25) we get 


B B. 
tan — = ——*—_ (32) 
ye Zs PRIA /v 
The leakage rate is thus 
= 5B, 
y=65,—-——_* (33) 


é 
nd the time constant 


SP Rat af ¢ epee ea (34) 
tan B,/Xs 
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This is a more exact formula than that given in eq. (1). As the ratio B, pS, is a 
small quantity in pulsed source experiments, we can make a series development 
and find 


Ap =U Let 


=~ 22 2, ¢ 4 
ae Be eae get (35) 


LB 15S Sib Ss 


This series corresponds to that obtained for Z when the absorption is not small 
(D. 5.2). Usually for a moderator only the first term in the L-series has to be 
taken into account. However, for a pulsed source measurement the term Bj is 
also important as the measured time constant is fitted to a parabola in Bj and 
the coefficient for the Bt term is usually interpreted as the diffusion cooling coeffi- 
cient according to eq. (3). The order of magnitude of the BS term also may be 
of interest in assessing the validity of the assumed parabolic dependence. 

Thus eq. (1) for the time constant should also contain a term d Bj and the 
coefficient for B‘ in eq. (3) is then C-d. This means that all the diffusion cooling 
coefficients measured up to now should be increased by d.1 Before calculating the 
magnitude of this correction we have to study the influence of anisotropic scattering. 


The relation between time constant and buckling for anisotropic scattering 


When the scattering is anisotropic it is convenient to express the scattering 
function in eq. (21) as an expansion in Legendre polynomials (D. 17.1) 


ties sts , } 
ae 2 (2 +1) bp Pn (Q + 2"), (36) 


{(Q'+Q)= 


where b,=1. For simplicity we study here only the case of plane geometry and 
the Boltzmann equation then takes the form 


Ow (x, u x; OR nd 

u vel ve Yo Yee (2 +1) bn Wn (x) Pn (ee), (37) 
1 ao 

where wy (@, u)= 7 > (204+ 1) yn (x) Pr (u) (38) 
n=0 


and yz is thé x-component of &. 

Just as for isotropic scattering a relation between-c, L and I can be derived 
(D. 17.3) but in the form of a determinantal equation. In the first order approxi- 
mation it is found that if the transport cross section defined by 


Dir = Ds (1-2) (39) 


is introduced instead of the scattering cross section, the same expression for the 
diffusion length can be used. We can thus expect the expression for the time 
constant to be correct in the B? term if the transport cross section is introduced, 

1 Kazarnovsky, Stepanoy and Shapiro have mentioned this correction for the case of the P; 


approximation in plane geometry (Second Conf. on the Peaceful Uses of Atomic Energy, Paper 
P/2148, Geneva 1958). 


154 


ARKIV FOR FYSIK. Bd 15 nr 12 


However, due to the importance of the Bj term for the interpretation of the 
diffusion cooling effect a higher order approximation is needed. 

As the calculation and solution of the determinantal equation is rather labo- 
rious, the spherical harmonics method was used instead. Multiplication of eq. (37) 
with (2n+1)P,(u) and integration over all w gives the following system of diffe- 
rential equations 


1—cb, 
l 


(n +1) phir (a)+nph_1 (a) + (2041) yn()=0 (n=0,1,2...). (40) 


We seek solutions of the form 


Wn (2) =9,e"* (41) 


and find that the following, simpler determinant must vanish: 


l-—c tal 0 0 
vo 3 (1603) Zial 0 
0 2400-6 (V—cb,) 3ial he , 
0 0 os erat en eee ee cee) 


Substituting /* and c* from eqs. (24) and (25) and identifying the positive value 

for x with Bj gives the following result in the P, approximation 
6 Bs (,, Bo_1+4b,—5b, 

3 Der 15 Dt i! a : 


Ay =U Xiet (43) 


The P,; approximation leads to the same expression if b, and 6, are assumed neg- 
ligible. : 

When the spherical harmonics method is applied to a spherically symmetrical 
system the relation between the time constant and the buckling is given by the 
same determinant (D. 11.1). Therefore eq. (43) can be regarded as a good approxi- 
mation for the geometries used in pulsed source experiments. 

As the Bj term turns out to be very small in eq. (35) compared to the uncer- 
tainties in the dependence of the diffusion constant on the buckling it is not ne- 
cessary to calculate the influence of the anisotropy on that term. 


The corrections to the diffusion cooling coefficient 


We can now estimate corrections to be applied to diffusion cooling coefficients 
which have been calculated from measurements by use of eq. (3). The corrections 
are given by the last term in eq. (43) which, we may recall, was derived under 
the assumption that the angular distribution in the scattering process is inde- 
pendent of neutron energy. 

It is well known from theory and experiment that this assumption is not valid 
for most moderators, so the coefficients 6, and b, have to be interpreted as some 
average values over the neutron spectrum. In lack of detailed information we 
ssume that b,=0, and the correction d, which is always positive, is then cal- 
ulated from the formula 
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Drie 
a= 0 “tr 
15 


(1+ 46,). (44) 


The averages used for b, of H,O, CH, and D,O have been obtained from eq. (39) 


where ©, has been calculated from measured curves and &;, from measured values 
of l,, taken from the literature. For H,O and D,O the curves given in BNL-325 
were used, but no curve was available for CH,. Therefore the cross section of 
polythene was measured using the slow and fast choppers at the Swedish Rl 
reactor.! For the moderators Be, BeO and graphite we made the reasonable as- 
sumption 6, =0. 

The calculated corrections, d, are listed in table 1 together with adopted values 
for the parameters used and for the coefficient C—d of the Bs term in eq. (3). 
The errors given do not include the possible uncertainty from higher terms in 
eq. (2). As so many approximations have been made in the evaluation the result 
must be considered only to give an estimate of the magnitude, especially for the 
three first moderators. However, it is obvious that the corrections are by no means 
negligible. 


Table 1. 
te ee SS 
U d C-d References for 
0 
Medges tons S106 cm?/see em | b | 10° cm4/sec | 10° cm*/sec C-d 
| 

H,O 3.52 0.425 0.25 0.0085 | 0.036+0.015| [1, 2, 4, 5, 8] 
CH, 2.66 0.323 0.25 0.0037 0.012 + 0.005 | [8, 9] 
D,O 20.8 2.51 0.12 1.3 3.54 0.8 [10] 
Be 12.2 1.48 0 0.18 341 [2, 5] 
BeO 11.8 1.43 0 0.16 4+] (3, 7] 
Graphite 21.2 2.56 0 0.93 14+2 [2, 6] 


The boundary conditions 


In the preceding sections we have obtained relations between the time constant 
and the buckling. The next problem is to investigate how the geometric buckling 
is calculated from the dimensions of a body. This usually involves the introduction 
of the eatrapolated endpoint x) which is defined as the distance beyond the surface 
of the body.at which the asymptotic neutron flux vanishes (D. 6.2). This is not 
the same as the linear extrapolation length, which is the ratio between the asymp- 
totic neutron flux and its derivative at the surface, and is only equal to a when 
c=1. 

Because of the small dimensions of the moderators used in pulsed neutron 
source experiments great care must be taken to find the right value of ay. 
large error here will cause a large error in the buckling of small systems and 
hence also in the result obtained for the diffusion cooling coefficient. For the 
stationary case formulae for a) have been given for anisotropic scattering and for 
various values of c (D. 17.4) but they are not directly applicable to the small 
geometries used in pulsed source experiments. 


1 My thanks are due to Dr. K.-E. Larsson and Mr. EK. Johansron for arranging these 
measurements. 
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Therefore, a calculation was made for an infinite slab of thickness 2a using 
the P, approximation of the spherical harmonics method. We start then from 
the system of equations (40) assuming the coefficients b, and above to be zero. 
For c and / the expressions (24) and (25) are used. With the origin in the centre 
of the slab the solution for wy(x) is found to be of the form 


Wo (x) = A, cos ax+ A, cosh Bx, (45) 


where using eq. (42) for c* close to unity « and f can be expressed as 


a21*2 = 3 (c* — 1) (1—c*},) (1+ $ (c*—1)), is) 
sores (a=) 


From the system of equations (40) similar expressions for ,, y, and w, are easily 
obtained. 
The most appropriate boundary conditions are those named after Marshak 


(D. 10.3) 


0 0 
Iv bw) wdu=O, lv bw) wedu=0. (48) 
=1 =i 


Introducing the series expansion of eq. (38) they can also be expressed as 


ae ea Mote av (a)=0, (49) 
3 5 2 
Pot) =, (a) +2 Y2 (a) ~= Ye (a) =0. (50) 


Using these equations and approximating, assuming c* to be close to unity, an 
expression for tan aa is derived. 

The asymptotic neutron flux is represented by the cosine term in eq. (45). The 
extrapolated endpoint is therefore defined by 


a (a+ %) => (51) 


and x is just the geometric buckling Bj. In the first approximation we find 
a = 0.7051 [1 — B31z. (0.0256 — 0.2825 b, + 0.0893 b,2)]. (52) 
tr 


The factor 0.7051 is close to the correct value 0.7104, which shows that the cal- 
culation is rather accurate. 


For 6,=0, —%/ltr = 0.7051 (1 — 0.0256 B5 1:2) 
For b,=0.25, x9/liz = 0.7051 (1 + 0.0394 B5 1,). 


_ These figures show that for systems used in pulsed source experiments, where 
Bi lz. is usually smaller than 0.1, the resulting correction to the extrapolated end- 
int is so small that it can be neglected. The same conclusion was also reached 
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in a similar calculation for a sphere. This is in agreement with the numerical 
results mentioned in reference 1, which were obtained with the assumption of 
isotropic scattering. 

The finite geometry may cause indirectly another uncertainty in the choice 
of a suitable extrapolated endpoint. The diffusion cooling effect causes a change 
in the neutron spectrum which as a first approximation can be described as a 
decrease in the temperature of the Maxwell spectrum of the neutrons [15]. The 
average scattering cross section defined in eq. (23) and hence the transport mean 
free path may therefore vary with buckling. For water the decrease in neutron 
temperature is of the order of 70° per unit buckling (measured in em ”) [1] [14]. 
As the transport mean free path is rather accurately proportional to the square 
root of the neutron temperature the extrapolated endpoint will be about 10 per 
cent smaller at a buckling of 0.8 em~”. Corrections are therefore necessary. 


Conclusions 


Even if some results have been obtained for the problems discussed here no 
full understanding has been reached. More experimental data on the scattering 
functions of moderators are needed and also further theoretical investigations. 
Especially, the possibility of a spatial dependence of the neutron velocity distri- 
bution should be studied, as such an effect might be of importance in the small 
systems treated here. 
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